package org.jmap.core;

import java.util.HashMap;
import java.util.Map;

import org.osgeo.proj4j.units.Unit;
import org.osgeo.proj4j.units.Units;

public class Util {
	
	public static long lastSequentialId = 0;

	private static Map<String, Double> INCHES_PER_UNIT = new HashMap<String, Double>() {{
		put("inches", 1.0);
		put("ft", 12.0);
		put("mi", 63360.0);
		put("m", 39.3701);
		put("km", 39370.1);
		put("dd", 4374754d);
		put("yd", 36d);
		put("degrees", get("dd"));
		put("in", get("inches"));
		put("nmi", get("m") * 1852);
	}};
	
	public static Map<String, Double> VINCENTY_CONSTANTS = new HashMap<String, Double>() {{
		    put("a", 6378137d);
		    put("b", 6356752.3142d);
		    put("f", 1/298.257223563d);
	}};
	
	public static final double METERS_PER_INCH = 0.02540005080010160020;
	public static final double DOTS_PER_INCH = 72.0;
	
	public static double normalizeScale(double scale) {
		if (scale > 1.0) {
			return 1.0/scale;
		} else {
			return scale;
		}
	}
	
	public static String generateUniqueId(String prefix) {
	    if (prefix == null) {
	        prefix = "id_";
	    }
	    
	    lastSequentialId += 1; 
	    return prefix + lastSequentialId;
	}
	
	/**
	 * gets the resolution for a given scale
	 * 
	 * @param scale the scale to find resolution for
	 * @param units the units to use, DEGREES by default
	 * @return the resolution at the given scale
	 */
	public static double getResolutionFromScale(double scale, Unit units) {
	    double resolution;
	    
        if (units == null) {
            units = Units.DEGREES;
        }
        
        double normScale = normalizeScale(scale);
        resolution = 1 / (normScale * INCHES_PER_UNIT.get(units.plural) * DOTS_PER_INCH);        
	    
	    return resolution;
	}
	
	/**
	 * gets the scale for a given resolution
	 * 
	 * @param the resolution to find the scale for
	 * @param units Index into INCHES_PER_UNIT hashmap.
	 *                  Default is degrees
	 * 
	 * @Returns
	 * The corresponding scale given passed-in resolution and unit 
	 *         parameters.
	 */
	public static double getScaleForResolution(double resolution, Unit units) {
	    if (units == null) {
	        units = Units.DEGREES;
	    }

	    double scale = resolution * INCHES_PER_UNIT.get(units.plural) * DOTS_PER_INCH;
	    return scale;
	}
	
	public static double degreesToRadians(double x) {
		return x * Math.PI/180;
	}
	
	
	public static double geodeticDistance(MapPoint p1, MapPoint p2) {
		double a = VINCENTY_CONSTANTS.get("a");
		double b = VINCENTY_CONSTANTS.get("b");
		double f = VINCENTY_CONSTANTS.get("f");

	    double L = degreesToRadians(p2.lon - p1.lon);
	    double U1 = Math.atan((1-f) * Math.tan(degreesToRadians(p1.lat)));
	    double U2 = Math.atan((1-f) * Math.tan(degreesToRadians(p2.lat)));
	    double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
	    double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
	    double lambda = L, lambdaP = 2*Math.PI;
	    double iterLimit = 20d;
	    double sigma = 0, cosSqAlpha = 0, sinLambda = 0, sinSigma = 0, cosLambda = 0, cosSigma = 0, cos2SigmaM = 0;
	    
	    while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0) {
	        sinLambda = Math.sin(lambda); 
	        cosLambda = Math.cos(lambda);
	        sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
	        (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
	        if (sinSigma==0) {
	            return 0;  // co-incident points
	        }
	        cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
	        sigma = Math.atan2(sinSigma, cosSigma);
	        double alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
	        cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
	        cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
	        double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
	        lambdaP = lambda;
	        lambda = L + (1-C) * f * Math.sin(alpha) *
	        (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
	    }
	    
	    if (iterLimit==0) {
	        return Double.NaN;  // formula failed to converge
	    }
	    double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
	    double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
	    double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
	    double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
	        B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
	    double s = b*A*(sigma-deltaSigma);
	    double d = s/1000; // round to 1mm precision
	    return d;
	}
	
}
